{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to atomman: Free surface generator\n",
    "\n",
    "__Lucas M. Hale__, [lucas.hale@nist.gov](mailto:lucas.hale@nist.gov?Subject=ipr-demo), _Materials Science and Engineering Division, NIST_.\n",
    "    \n",
    "[Disclaimers](http://www.nist.gov/public_affairs/disclaimer.cfm) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Introduction<a id='section1'></a>\n",
    "\n",
    "This Notebook outlines the FreeSurface class that allows for the easy generation of planar-sliced free surface atomic configurations based on a unit cell and Miller (hkl) or Miller-Bravais (hkil) crystal plane indices for any crystal unit cell.\n",
    "\n",
    "*New as of version 1.3.2*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Library imports**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "atomman version = 1.4.0\n",
      "Notebook executed on 2021-08-05\n"
     ]
    }
   ],
   "source": [
    "# Standard Python libraries\n",
    "import datetime\n",
    "\n",
    "# http://www.numpy.org/\n",
    "import numpy as np\n",
    "\n",
    "# https://github.com/usnistgov/atomman\n",
    "import atomman as am\n",
    "import atomman.unitconvert as uc\n",
    "\n",
    "# Show atomman version\n",
    "print('atomman version =', am.__version__)\n",
    "\n",
    "# Show date of Notebook execution\n",
    "print('Notebook executed on', datetime.date.today())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Theory<a id='section2'></a>\n",
    "\n",
    "A planar-sliced free surface atomic system can be constructed by\n",
    "\n",
    "1. Generating a bulk crystalline system with all periodic boundaries.\n",
    "2. Slicing the system along a geometrical plane creating two ideal planar-sliced free surfaces.\n",
    "\n",
    "While these steps are conceptually simple, it is still tricky to create a routine that can successfully work for any crystal plane in any crystal. The main challenge in developing a generalized free surface routine is working with the multiple different axes systems that describe the defect and the atomic configurations.\n",
    "\n",
    "Planar free surfaces are typically defined as crystallographic $(hkl)$ Miller or $(hkil)$ Miller-Bravais planes relative to conventional unit cells with lattice box vectors of $a_{unit}$, $b_{unit}$, and $c_{unit}$. [Sun and Cedar](https://doi.org/10.1016/j.susc.2013.05.016) showed an algorithm that can be easily used to identify the Cartesian slip plane normal for any $(hkl)$ plane in any crystal. In atomman, this is handled by atomman.tools.miller.plane_crystal_to_cartesian.\n",
    "\n",
    "Using the Cartesian plane normal, a search is then performed to identify three $[uvw]$ crystal vectors corresponding to\n",
    "\n",
    "- The crystal vector with absolute $h$, $k$, and $l$ values less than some search maximum that is closest to the plane normal.\n",
    "- A crystal vector in the slip plane that has the shortest length.\n",
    "- A second crystal vector in the slip plane that has the shortest length without being parallel to the previous vector. \n",
    "\n",
    "Uniqueness is supported by only picking vectors that form a right-handed system and which have the smallest angle between the two in-plane vectors.\n",
    "\n",
    "A rotated cell is then constructed in which the rotated cell's box vectors $a_{rot}$, $b_{rot}$, and $c_{rot}$ correspond to the three $[uvw]$ Miller vectors of the unit cell. The specific alignment of the rotated cell depends on which box vector the free surface cut will be made along, i.e. the \"cut box vector\"\n",
    "\n",
    "- The cut box vector is aligned with the out-of-plane crystal vector.\n",
    "- The Cartesian system for the rotated cell is defined such that the cut box vector is the only box vector with a component along one of the three x, y, or z Cartesian axes.\n",
    "\n",
    "This choice means that the free surface can be inserted into the system easily by making the boundary conditions along the cut box vector non-periodic.  As the other two box vectors are in-plane vectors, the resulting free surface plane is the plane of interest. Also, the choice of Cartesian system mentioned will result in the plane normal coinciding with one of the three Cartesian axes.\n",
    "\n",
    "__Note__: Rotations in atomman create systems that are compatible with the box vector limitations imposed by LAMMPS\n",
    "\n",
    "- $a = [lx, 0, 0]$\n",
    "- $b = [xy, ly, 0]$\n",
    "- $c = [xz, yz, lz]$\n",
    "\n",
    "As such, $c$ is the preferred cut box vector as it is the only one with a component in the z direction, and therefore will work for any crystal system.\n",
    "\n",
    "The resulting free surface system can be constructed from the rotated cell by \n",
    "\n",
    "1. Creating a supercell by replicating along the rotated cell box vectors\n",
    "2. Rigidly shifting all atoms in the system so that the box edge where the cut will occur falls between two atomic planes\n",
    "3. Creating the free surface by either making the cut box vector boundary condition non-periodic, or by increasing the cut box vector to effectively insert a region of vacuum.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Class basics<a id='section3'></a>\n",
    "\n",
    "The class is set up to generate a free surface system in two steps:\n",
    "\n",
    "1. Initialize the class and specify the (hkl) surface and unit cell structure to use. The proper rotated cell will be generated as well as a list of shifts that can be applied.\n",
    "2. Build the free surface system by selecting the shift value corresponding to the termination planes of interest and the size of the final system.\n",
    "\n",
    "### 3.1. FreeSurface initialization\n",
    "\n",
    "- __hkl__ (*array-like object*) The free surface plane to generate expressed in either 3 indices Miller (hkl) format or 4 indices Miller-Bravais (hkil) format.\n",
    "- __ucell__ (*atomman.System*) The unit cell to use in generating the system.\n",
    "- __cutboxvector__ (*str, optional*) Specifies which of the three box vectors corresponds to the out-of-plane vector.  Default value is c.\n",
    "- __maxindex__ (*int, optional*) Max uvw index value to use in identifying the best uvw set for the out-of-plane vector.  If not given, will use the largest absolute index between the given hkl and the initial in-plane vector guesses.\n",
    "- __conventional_setting__ (*str, optional*) Allows for rotations of a primitive unit cell to be determined from (hkl) indices specified relative to a conventional unit cell.  Allowed settings: 'p' for primitive (no conversion), 'f' for face-centered, 'i' for body-centered, and 'a', 'b', or 'c' for side-centered.  Default behavior is to perform no conversion, i.e. take (hkl) relative to the given ucell.\n",
    "- __tol__ (*float, optional*) Tolerance parameter used to round off near-zero values.  Default value is 1e-8.\n",
    "\n",
    "\n",
    "### 3.2. FreeSurface.surface()\n",
    "\n",
    "- __shift__ (*array-like object, optional*) Applies a shift to all atoms. Different values allow for free surfaces with different termination planes to be selected.\n",
    "- __vacuumwidth__ (*float, optional*) If given, the free surface is created by modifying the system's box to insert a region of vacuum with this width. This is typically used for DFT calculations where it is computationally preferable to insert a vacuum region and keep all dimensions periodic.\n",
    "- __sizemults__ (*list or tuple, optional*) The three System.supersize multipliers [a_mult, b_mult, c_mult] to use on the rotated cell to build the final system. Note that the cutboxvector sizemult must be an integer and not a tuple.  Default value is [1, 1, 1].\n",
    "- __minwidth__ (*float, optional*) If given, the sizemult along the cutboxvector will be selected such that the width of the resulting final system in that direction will be at least this value. If both sizemults and minwidth are given, then the larger of the two in the cutboxvector direction will be used. \n",
    "- __even__ (*bool, optional*) A True value means that the sizemult for cutboxvector will be made an even number by adding 1 if it is odd.  Default value is False."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Examples<a id='section4'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1. Conventional bcc unit cell\n",
    "\n",
    "#### 4.1.1. Define a bcc conventional cell and initialize a FreeSurface"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = uc.set_in_units(2.866, 'angstrom')\n",
    "box = am.Box.cubic(a)\n",
    "\n",
    "atoms = am.Atoms(pos=[[0.0, 0.0, 0.0], [0.5, 0.5, 0.5]])\n",
    "\n",
    "ucell = am.System(atoms=atoms, box=box, scale=True, symbols='Fe')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "plane_hkl = [1, 1, 0]\n",
    "\n",
    "fs = am.defect.FreeSurface(plane_hkl, ucell)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 4.1.2. Check the identified rotated system\n",
    "\n",
    "The FreeSurface class has a number of attributes that can be used to check the solution.\n",
    "\n",
    "- __uvws__ are the three conventional Miller vectors used as the box vectors of the rotated system.  \n",
    "- __rcell__ is the rotated cell.\n",
    "- __transform__ is the Cartesian transformation matrix associated with going from ucell to rcell.\n",
    "- __rcellwidth__ is the width of the rotated cell in the direction normal to the hkl plane."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.  0.  1.]\n",
      " [ 1. -1.  0.]\n",
      " [ 1.  1.  0.]]\n"
     ]
    }
   ],
   "source": [
    "print(fs.uvws)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "avect =  [ 2.866,  0.000,  0.000]\n",
      "bvect =  [ 0.000,  4.053,  0.000]\n",
      "cvect =  [ 0.000,  0.000,  4.053]\n",
      "origin = [ 0.000,  0.000,  0.000]\n",
      "natoms = 4\n",
      "natypes = 1\n",
      "symbols = ('Fe',)\n",
      "pbc = [ True  True  True]\n",
      "per-atom properties = ['atype', 'pos']\n",
      "     id |   atype |  pos[0] |  pos[1] |  pos[2]\n",
      "      0 |       1 |   1.433 |   2.027 |   0.000\n",
      "      1 |       1 |   0.000 |   0.000 |   0.000\n",
      "      2 |       1 |   1.433 |   0.000 |   2.027\n",
      "      3 |       1 |   0.000 |   2.027 |   2.027\n"
     ]
    }
   ],
   "source": [
    "print(fs.rcell)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.00000000e+00  4.46536745e-18  1.00000000e+00]\n",
      " [ 7.07106781e-01 -7.07106781e-01  3.19640147e-16]\n",
      " [ 7.07106781e-01  7.07106781e-01 -9.73839539e-18]]\n"
     ]
    }
   ],
   "source": [
    "print(fs.transform)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.05313606976129\n"
     ]
    }
   ],
   "source": [
    "print(fs.rcellwidth)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 4.1.3. Build the free surface configuration\n",
    "\n",
    "Building the final free surface configuration requires specifying an atomic shift to position the fault plane between two atomic planes.  The FreeSurface.shifts attribute lists all identified shift positions for rcell that will place the geometrical plane halfway between atomic planes. A shift value is not automatically used as a default as different shifts can result in distinctly different termination planes being created. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.         0.         1.01328403]\n",
      " [0.         0.         3.03985207]]\n"
     ]
    }
   ],
   "source": [
    "print(fs.shifts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a given shift, a free surface configuration can then be retrieved using the surface method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "avect =  [ 2.866,  0.000,  0.000]\n",
      "bvect =  [ 0.000,  4.053,  0.000]\n",
      "cvect =  [ 0.000,  0.000, 12.159]\n",
      "origin = [ 0.000,  0.000,  0.000]\n",
      "natoms = 12\n",
      "natypes = 1\n",
      "symbols = ('Fe',)\n",
      "pbc = [ True  True False]\n",
      "per-atom properties = ['atype', 'pos']\n",
      "     id |   atype |  pos[0] |  pos[1] |  pos[2]\n",
      "      0 |       1 |   1.433 |   2.027 |   1.013\n",
      "      1 |       1 |   0.000 |   0.000 |   1.013\n",
      "      2 |       1 |   1.433 |   0.000 |   3.040\n",
      "      3 |       1 |   0.000 |   2.027 |   3.040\n",
      "      4 |       1 |   1.433 |   2.027 |   5.066\n",
      "      5 |       1 |   0.000 |   0.000 |   5.066\n",
      "      6 |       1 |   1.433 |   0.000 |   7.093\n",
      "      7 |       1 |   0.000 |   2.027 |   7.093\n",
      "      8 |       1 |   1.433 |   2.027 |   9.120\n",
      "      9 |       1 |   0.000 |   0.000 |   9.120\n",
      "     10 |       1 |   1.433 |   0.000 |  11.146\n",
      "     11 |       1 |   0.000 |   2.027 |  11.146\n"
     ]
    }
   ],
   "source": [
    "print(fs.surface(shift=fs.shifts[0], minwidth=10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After a free surface configuration has been created, the FreeSurface object will also have a __surfacearea__ attribute which gives the surface area for one of the two free surfaces created.\n",
    "\n",
    "Note: For free surface formation energy calculations you must divide by 2\\*surfacearea as two free surfaces are created."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "11.616287975935858\n"
     ]
    }
   ],
   "source": [
    "print(fs.surfacearea)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.2. Primitive bcc unit cell\n",
    "\n",
    "The algorithm relies on identifying the smallest in-plane Miller vectors related to the unit cell given. If the conventional unit cell is not primitive, then the identified vectors may not be the smallest lattice vectors possible, and therefore the system may be larger than needed.  Alternately, the primitive unit cell can be used to construct the free surface configuration by specifying the conventional_setting parameter\n",
    "\n",
    "- Give ucell as the primitive unit cell.\n",
    "- Use conventional_setting to indicate the Bravais lattie centering: 'p' for primitive (no conversion), 'f' for face-centered, 'i' for body-centered, and 'a', 'b', or 'c' for side-centered.\n",
    "- Then, (hkl) can be specified relative to the conventional unit cell, but used to operate on the primitive unit cell."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 4.2.1. Define a primitive bcc unit cell and initialize a FreeSurface"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "box = am.Box.trigonal(a = a * 3**0.5 / 2, alpha = 109.47122063449069)\n",
    "\n",
    "atoms = am.Atoms(pos=[[0.0, 0.0, 0.0]])\n",
    "\n",
    "ucell = am.System(atoms=atoms, box=box, scale=True, symbols='Fe')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note we are giving the same plane_hkl value as above (relative to the conventional cubic cell), but specifying the conventional setting. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "plane_hkl = [1, 1, 0]\n",
    "setting = 'i'\n",
    "fs = am.defect.FreeSurface(plane_hkl, ucell, conventional_setting=setting)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 4.1.2. Check the identified rotated system\n",
    "\n",
    "Using the primitive unit cell allows for the algorithm to find shorter lattice vectors than can be identified with the conventional unit cell. The identified uvws (relative to the conventional cell) may then be fractions instead of integers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-0.5  0.5 -0.5]\n",
      " [-0.5  0.5  0.5]\n",
      " [ 1.   1.   0. ]]\n"
     ]
    }
   ],
   "source": [
    "print(fs.uvws)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For rcell, the system is half the size but still has the same width."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "avect =  [ 2.482,  0.000,  0.000]\n",
      "bvect =  [ 0.827,  2.340,  0.000]\n",
      "cvect =  [ 0.000,  0.000,  4.053]\n",
      "origin = [ 0.000,  0.000,  0.000]\n",
      "natoms = 2\n",
      "natypes = 1\n",
      "symbols = ('Fe',)\n",
      "pbc = [ True  True  True]\n",
      "per-atom properties = ['atype', 'pos']\n",
      "     id |   atype |  pos[0] |  pos[1] |  pos[2]\n",
      "      0 |       1 |   0.000 |   0.000 |   0.000\n",
      "      1 |       1 |   1.655 |   1.170 |   2.027\n"
     ]
    }
   ],
   "source": [
    "print(fs.rcell)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 4.1.3. Build the free surface configuration\n",
    "\n",
    "The resulting system is then more friendly to DFT by having fewer atoms in the periodic directions. Here, we also use the vacuumwidth parameter to insert a region of vacuum as that is also computationally favorable for DFT calculations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "avect =  [ 2.482,  0.000,  0.000]\n",
      "bvect =  [ 0.827,  2.340,  0.000]\n",
      "cvect =  [ 0.000,  0.000, 27.159]\n",
      "origin = [ 0.000,  0.000, -7.500]\n",
      "natoms = 6\n",
      "natypes = 1\n",
      "symbols = ('Fe',)\n",
      "pbc = [ True  True False]\n",
      "per-atom properties = ['atype', 'pos']\n",
      "     id |   atype |  pos[0] |  pos[1] |  pos[2]\n",
      "      0 |       1 |   0.000 |   0.000 |   1.013\n",
      "      1 |       1 |   1.655 |   1.170 |   3.040\n",
      "      2 |       1 |   0.000 |   0.000 |   5.066\n",
      "      3 |       1 |   1.655 |   1.170 |   7.093\n",
      "      4 |       1 |   0.000 |   0.000 |   9.120\n",
      "      5 |       1 |   1.655 |   1.170 |  11.146\n"
     ]
    }
   ],
   "source": [
    "print(fs.surface(shift=fs.shifts[0], minwidth=10, vacuumwidth=15))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5.808143987967931\n"
     ]
    }
   ],
   "source": [
    "print(fs.surfacearea)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.3. hcp unit cell\n",
    "\n",
    "Hexagonal systems can also be worked with using Miller-Bravais planes and vectors."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 4.3.1. Define an hcp unit cell and initialize a FreeSurface"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "avect =  [ 3.185,  0.000,  0.000]\n",
      "bvect =  [-1.593,  2.759,  0.000]\n",
      "cvect =  [ 0.000,  0.000,  5.185]\n",
      "origin = [ 0.000,  0.000,  0.000]\n",
      "natoms = 2\n",
      "natypes = 1\n",
      "symbols = ('Mg',)\n",
      "pbc = [ True  True  True]\n",
      "per-atom properties = ['atype', 'pos']\n",
      "     id |   atype |  pos[0] |  pos[1] |  pos[2]\n",
      "      0 |       1 |   0.000 |   0.000 |   0.000\n",
      "      1 |       1 |   0.000 |   1.839 |   2.593\n"
     ]
    }
   ],
   "source": [
    "a = uc.set_in_units(3.1853, 'angstrom')\n",
    "c = uc.set_in_units(5.1853, 'angstrom')\n",
    "box = am.Box.hexagonal(a, c)\n",
    "\n",
    "atoms = am.Atoms(pos = [[0.0, 0.0, 0.0], [1/3, 2/3, 0.5]])\n",
    "\n",
    "ucell = am.System(atoms=atoms, box=box, scale=True, symbols='Mg')\n",
    "print(ucell)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "plane_hkl = [1, 0, -1, 0]\n",
    "fs = am.defect.FreeSurface(plane_hkl, ucell)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 4.1.2. Check the identified rotated system\n",
    "\n",
    "If the plane is defined using Miller-Bravais indices, then the uvws will be converted to them as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-0.33333333  0.66666667 -0.33333333  0.        ]\n",
      " [ 0.          0.         -0.          1.        ]\n",
      " [ 0.66666667 -0.33333333 -0.33333333  0.        ]]\n"
     ]
    }
   ],
   "source": [
    "print(fs.uvws)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "avect =  [ 3.185,  0.000,  0.000]\n",
      "bvect =  [ 0.000,  5.185,  0.000]\n",
      "cvect =  [-1.593,  0.000,  2.759]\n",
      "origin = [ 0.000,  0.000,  0.000]\n",
      "natoms = 2\n",
      "natypes = 1\n",
      "symbols = ('Mg',)\n",
      "pbc = [ True  True  True]\n",
      "per-atom properties = ['atype', 'pos']\n",
      "     id |   atype |  pos[0] |  pos[1] |  pos[2]\n",
      "      0 |       1 |   1.593 |   0.000 |   2.759\n",
      "      1 |       1 |   1.593 |   2.593 |   0.920\n"
     ]
    }
   ],
   "source": [
    "print(fs.rcell)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.         0.         0.91951692]\n",
      " [0.         0.         2.29879228]]\n"
     ]
    }
   ],
   "source": [
    "print(fs.shifts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 4.1.3. Build the free surface configuration\n",
    "\n",
    "If you give both sizemults and minwidth, then the larger of the two will be used for the determining the final system width."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "avect =  [ 3.185,  0.000,  0.000]\n",
      "bvect =  [ 0.000,  5.185,  0.000]\n",
      "cvect =  [-9.556,  0.000, 16.551]\n",
      "origin = [ 0.000,  0.000,  0.000]\n",
      "natoms = 12\n",
      "natypes = 1\n",
      "symbols = ('Mg',)\n",
      "pbc = [ True  True False]\n",
      "per-atom properties = ['atype', 'pos']\n",
      "     id |   atype |  pos[0] |  pos[1] |  pos[2]\n",
      "      0 |       1 |  -1.593 |   0.000 |   3.678\n",
      "      1 |       1 |   1.593 |   2.593 |   1.839\n",
      "      2 |       1 |  -3.185 |   0.000 |   6.437\n",
      "      3 |       1 |   0.000 |   2.593 |   4.598\n",
      "      4 |       1 |  -4.778 |   0.000 |   9.195\n",
      "      5 |       1 |  -1.593 |   2.593 |   7.356\n",
      "      6 |       1 |  -6.371 |   0.000 |  11.954\n",
      "      7 |       1 |  -3.185 |   2.593 |  10.115\n",
      "      8 |       1 |  -7.963 |   0.000 |  14.712\n",
      "      9 |       1 |  -4.778 |   2.593 |  12.873\n",
      "     10 |       1 |   0.000 |   0.000 |   0.920\n",
      "     11 |       1 |  -6.371 |   2.593 |  15.632\n"
     ]
    }
   ],
   "source": [
    "print(fs.surface(shift=fs.shifts[0], sizemults=(1, 1, 6), minwidth=10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also, note that the two shifts result in different termination planes. Both can be retrived from the same FreeSurface object by simply giving surface() a different shift value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "avect =  [ 3.185,  0.000,  0.000]\n",
      "bvect =  [ 0.000,  5.185,  0.000]\n",
      "cvect =  [-9.556,  0.000, 16.551]\n",
      "origin = [ 0.000,  0.000,  0.000]\n",
      "natoms = 12\n",
      "natypes = 1\n",
      "symbols = ('Mg',)\n",
      "pbc = [ True  True False]\n",
      "per-atom properties = ['atype', 'pos']\n",
      "     id |   atype |  pos[0] |  pos[1] |  pos[2]\n",
      "      0 |       1 |  -1.593 |   0.000 |   5.057\n",
      "      1 |       1 |  -1.593 |   2.593 |   3.218\n",
      "      2 |       1 |  -3.185 |   0.000 |   7.816\n",
      "      3 |       1 |  -3.185 |   2.593 |   5.977\n",
      "      4 |       1 |  -4.778 |   0.000 |  10.574\n",
      "      5 |       1 |  -4.778 |   2.593 |   8.735\n",
      "      6 |       1 |  -6.371 |   0.000 |  13.333\n",
      "      7 |       1 |  -6.371 |   2.593 |  11.494\n",
      "      8 |       1 |  -7.963 |   0.000 |  16.092\n",
      "      9 |       1 |  -7.963 |   2.593 |  14.253\n",
      "     10 |       1 |   0.000 |   0.000 |   2.299\n",
      "     11 |       1 |   0.000 |   2.593 |   0.460\n"
     ]
    }
   ],
   "source": [
    "print(fs.surface(shift=fs.shifts[1], sizemults=(1, 1, 6), minwidth=10))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
